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SOLUTION  OF  THE  ABEL  INTEGRAL  TRANSFORM  FOR  A 

CYLINDRICAL  LUMINOUS  REGION  WITH  OPTICAL 

DISTORTIONS  AT  ITS  BOUNDARY* 

Earl  R.    Mosburg,    Jr.    and  Matthew  S.    Lojko 

The  use  of  orthogonal  polynomial  expansions  in 
the  calculation  of  the  Abel  integral  transform  is  dis- 
cussed.    Particular  attention  is  directed  to  the  effects 
of  optical  and  instrumental  distortions  when  the  lumi- 
nous region  is  contained  by  a  cylindrical  glass  tube. 
An  easily  calculable  solution  of  the  Abel  integral  is 
presented  which  reduces  the  effect  of  such  distortions 
by  employing  a  weighting  function  which  has  a  maxi- 
mum at  the  center  and  vanishes  at  the  boundary.      This 
approach  results  in  a  more  accurate  solution  of  the 
Abel  integral  transform  in  the  case  where  significant 
optical  and  instrumental  distortions  are  present  near 
the  boundary  of  the  luminous  region. 

Key  Words:     Abel  transform,    Abel  inversion,    plasma 
diagnostics,    emissivity  profile,    radiance 
profile. 

INTRODUCTION 
In  order  to  obtain  the  radial  distribution  of  volume  light 
emissivity  within  a  cylindrical,    non-absorbing  luminous  region,    we 
must  solve  the  Abel  integral  transform  using  the  projected  brightness 
profile,    measured  by  scanning  the  detector  in  a  direction  perpendicular 
to  the  axis  of  the  tube.      If  the  projected  brightness  profile  is  f(x),    where 
x  is  the  ratio  of  the  distance  off  axis  to  the    radius  of  the  luminous 
region,    and  if  g(r)  is  the  corresponding  volume  emissivity  distribution, 
where     r    is  the  normalized  radius,    then  the  Abel  integral  transforms 
can  be  written  as 

*Work  supported  in  part  by  the  Advanced  Research  Projects  Agency. 


/    \  2  d        pi        £(x)x  dx 

g(r)  =   -  —     "rdr"    Jr     yj     1      1  (la) 


x    -r 


and 


r,    »       r*  1    g(r)  r  dr  ,. ,  . 

£(x)  =  L  fi  \      •  (lb) 


'r    -x 


or  alternatively  in  the  forms 

2         pi      f'(x)  dx 


and 


(r)=   -~-  Jr  J    2''    2  (1C) 


x   -r 


/V 


r     -x     g'(r)dr  (Id) 


where  the  primes  indicate  differentiation.      A  direct  numerical  solution 
for  g(r)  using  measured  values  of  f(x)  in  Eq.    (la)  or  Eq.    (lc)  is  subject 
to  considerable  error  due  to  the  behavior  of  the  denominator  in  the 
integrand  and  to  the  necessity  for  numerical  differentiation.      These 
difficulties  are  considerably  alleviated  by  first  making  a  least  square 

fit  of    f(x)  to  a  power  series  expansion  as  described  by  Freeman  and 

1 

Katz. 

A  more  convenient  expansion  in  terms  of  orthogonal  polynomials 

2 
has  been  reported  by  Herlitz     using  Tchebycheff  polynomials  of  the 

3 
second  kind.     Popenoe  and  Shumaker     have  used  Herlitz' s  method  as 

well  as  an  expansion  in  terms  of  Legendre  polynomials.      The  use  of 

orthogonal  polynomial  expansions  is  equivalent  to  a  weighted  least 

squares  analysis.      But  here,    because  of  the  orthogonality  of  the  basis 

functions,    the  coefficients  can  be  independently  calculated.      Each 

coefficient  can  then  be  tested  for  statistical  significance  and  the 


expansion  appropriately  truncated  without  any  prior,    ad  hoc  decision 
about  the  number  of  terms  to  be  used.      The  weighting  function    w(x)/v(x) 
of  Eq.    (5)  is  determined  once  a  particular  series  of  orthogonal  poly- 
nomials is  chosen  for  the  expansion. 

Sufficient  attention  has  not,    however,    been  given  to  the  use  of 
orthogonal  polynomial  expansions  in  the  case  where  optical  or  instru- 
mental distortions  are  introduced  at  the  boundary  of  the  luminous 

region,    as  for  example,    by  the  presence  of  a  glass  container.      In  this 

4 
case,    distortions  due  to  the  scattering  and  uneven  refraction     of  light 

in  the  tube  walls  may  become  important.      These  effects  are  a  maximum 

near  x  =   1  where  a  near  grazing  angle  is  involved  in  the  measurement. 

Furthermore,    the  finite  size  of  the  spectrometer  slit  introduces  an 

averaging  over  the  normalized  spacial  resolution  function  of  the 

instrument,    R(x-£),    such  that  the  measured  curve  becomes  a  function, 

b.{Q),   where  ^ 

*C  +  D/2 

I  f(x)R(x-C)dx  (le) 

C  -  D/2 
and  ±D/2  are  the  limiting  values  of  (x-C)  for  which  there  is  appreciable 

contribution  to  the  integral.      The  projected  brightness  profile,    f(x), 

can,    in  principle,    be  recovered  from    h(£)  by  an  appropriate  inversion 

of  Eq.    (le),    but  residual  errors  will  be  present.      These  errors  will 

also  be  larger  near  x  =   1  where  the  differences  between  functions    h(£) 

and  f(x)  are  largest,    i.  e.  ,    where  the  second  derivative  of  h(£)  is  more 

important.      These  distortions  are  particularly  large  when  it  is  desired 

/       2 

to  invert  projected  profiles  approximating  f(x)  =Wl-x   ,    which  corresponds 

to  g(r)  =  constant.     Here  the  second  derivative  of  h(C)  is  even  larger 
near  the  boundary  and  the  luminosity  is  now  high  in  the  region  of 
maximum  distortion. 


We  wish  to  stress  at  this  point  that,    in  contrast  to  the  distortions, 
most  projected  brightness  profiles  of  experimental  interest  vanish  at 
x  =   1     and  exhibit  maxima  at  or  near  the  center  of  the  light  source.     It 
is  now  clear  that  in  order  to  reduce  the  effect  of  the  distortions,    we 
would  like  a  weighting  function  in  Eq.    (5)  which  vanishes  at  x  =   1     and 
exhibits  a  maximum  at  x  =   0. 

In  this  paper  we  restrict  our  choice  of  polynomial  to  the  general 
class  of  Ultraspherical    or  Gegenbauer  polynomials,    which  includes 
Legendre  and  Tchebycheff  polynomials  as  special  cases.      In  what 

follows  we  will  use  the  notation  of  the  Handbook  of  Mathematical 

5 
Functions.        We  expand    f(x)  in  terms  of  general  Gegenbauer  polynomials 

as 

N  .      (a) 

f(x)  =  E  a  v(x)C       (u(x))  (2) 

n  n 

n  =   0 

where  v(x)  is  some  shape  function  to  be  chosen  and  u(x)  is  some 
function  of  x.  Substituting  Eq.  (2)  into  Eq.  (la)  we  arrive  at  the 
expression 

7        N  v(x)C((X)(u(x))xdx 

g(r)=   -^-    E         a      -±-     J1    JL .     (3) 

tt  n     rdr       *r  /   c  c 

n=   0  Vx     -   r 

5 
When  the  orthonormalization  integral  for  the  Gegenbauer  polynomials 

is  written  in  the  form 


fb  w(x)  C(CC)  (u(x))  C(CC)  (u(x))  dx  =  h       6  ,  (4) 

Jan  m  na     nm 


then  multiplying  Eq.    (2)  by    — — r-     C         (u(x))  and  using  Eq.    (4), 

v(x)  m 

we  obtain 


1  pb            w(x)      (a).    .    ..                                              . 

a     =— —  f(x)          .     C       (u(x))  dx                                      (5) 

n        h  *»  a             v(x)        n 

na 

which  allows  the  calculation  of  the   coefficients  a     needed  in  Eq.    (3). 

n 

It  is  clear  that  the  function    f(x)  can  be  written  as  the  sum  of  an 

experimentally  significant  part,    f   (x),    and  a  part  due  to  optical  and 

instrumental  distortions,    f  n(x);  that  is, 

d 


f(x)=  f   (x)  +  f   (x).  (6) 

e  d 


In  many  cases  it  may  be  convenient  to  further  split  the  experimental 

part  into  an  easily  soluble  approximate  form,    f   (x),    and  a  relatively 

small  perturbation  to  this  form,    f   (x),    so  that 

P 

f(x)  =  f   (x)  +  [f    (x)  +  f    (x)]   =  f   (x)  +  f    (x)     .  (7) 

a  p  d  a  c 

Experimentally  the  terms  are  not,    of  course,    separable  from  prior 

knowledge,    but  we  may  arbitrarily  separate  out  the  approximate  function 

f   (x).      The  same  final  result  is  obtained  by  performing  the  Abel    inversion 

of  these  terms  separately  and  then  summing.      Note  that  the  polynomial 

expansion  used  for  f   (x)  need  not  be  the  same  as  that  used  for  the  two 

a 

other  terms     of  Eq.    (7). 

The  problem  then  becomes  one  of  reducing  the  effect  of  the 

distortion  contribution,    f  n(x),    in  treating  the  combined  contribution, 

d 

f(x)  of  Eq.    (6),    or  f   (x)  of  Eq.    (7).      When  choosing  a  specific  form  of 

c 

Eq.    (2)  we  wish  therefore  to  satisfy  two  requirements: 


(A)  Proper  weighting  factor.      The  weighting  factor  in  Eq.    (5) 

should  be  such  that  the  contribution  to  the  calculation  of  a     is  reduced 

n 

where  the  distortion  contribution,    f  n(x)  is  largest.      One  would  therefore 

d 

like  the  weighting  function  w(x)/v(x)  to  approach  zero  as  x-»l   and  be  a 
maximum  at  the  center. 

(B)  Ease  of  calculation.      In  principle,    Eq.    (3)  can  always  be 
evaluated  numerically,    but  the  number  of  integrals  that  must  be  cal- 
culated can  be  very  large.      To  illustrate  this  point,    if  we  wish  to 
calculate  g(r)  for  S  different  values  of  r,    then  the  number  of  integrals 
becomes  N(S+1)  where  N  is  the  number  of  terms  in  the  polynomial 
expansion.      For  convenience,    then,    Eq.    (3)  should  be  directly  integrable 
in  some  closed  form,    or,    this  failing,    it  should  be  easily  calculable 

as,  for  example,  by  a  recurrence  relation  between  the  integrals  of 
order  n+Z,  n+1,  and  n.  In  this  paper  we  have  settled  for  the  latter 
condition  in  order  to  satisfy  requirement    A    in  full. 

Once  the  form  of  Eq.    (2)  has  been  set,    the  weighting  factor  of 
requirement  A     and  the  integrability  or  non-integrability  of  Eq.    (3)  in 
closed  form  are  fully  determined.      Thus  the  simultaneous  satisfaction 
of  requirements     A  and  B    must  be  somewhat  fortuitous.      The  number  of 
solutions  of  the  Abel  integral  equation  in  terms  of  well-known  polynomials 
is  severely  limited.     In  rows  1  thru  4  of  Table  1,    we  show  the  weighting 
factors  and  solutions  of  the  integrals  of  Eq.    (3)  for  several  choices  of  the 
form  of  Eq.    (2)  where  the  function  v(x)  has  been  chosen  to  allow  evaluation 
of  these  integrals  in  closed  form.     All  of  these  solutions  can  be  derived 
by  proper  manipulation  of  the  general  equation 
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1  2  a "  2     (a )    2 

(1-x   )        C    (2x   -l)x  dx 


rdr     /  p I    ~~T 

7x   -r 


-^^fTfr(l-r    )         C  (2/-l)-C  (ZAl) 

r(a+l)  n  n-1 


We  are  not  aware  of  any  closed  form  solutions  which  are  not  specific 
cases  derivable  from  this  equation.      The  form  of  v(x)  in  these  solutions 
is  closely  related  to  the  normalization  function  w(x)  and  therefore  the 
weighting  function  w(x)/v(x)  cannot  be  arbitrarily  chosen.      None  of  the 
solutions  listed  has  a  weighting  function  of  the  form  we  desire.     In  this 
paper  we  present  a  solution  of  the  Abel  integral  equation  which  involves 
a  weighting  function  of  the  desired  type  and  consists  of  an  expansion  in 
terms  of  Gegenbauer  (a=2)  polynomials  (Table  1,  Row  5). 

THE  SOLUTION  USING  GEGENBAUER  (cc=2)  POLYNOMIALS 
If  we  choose  an  expansion  of  the  form 

f(x)=Za    J\    -    x2  c}2)  (x)  (9) 

n  2n 

then  Eq.    (5)becomes 

an=M2n+83)(2n+l)     T  ^  f(x)  (1   "  ^  C  Zn    <X>  dx  ■     ^ 

and  we  immediately  see  that  criterion    A     is     satisfied.      For  the 
sake  of  completeness,    the  first  few  polynomials  of  interest  are  given 
below. 

C(02)(x)  =   1  C^2)(x)  =  12  x  2-2  C^2)(x)  =  80  x  4-48  x  2  +3 

C,     (x)  =  448  x      -480  x  4+  120  x  2-4 


It   remains  to  evaluate  the  equation 


2    N 
g(r)=  -    —  Z 


_J_   _d_    o  1 

■*       A-~     .1 


1   -  x2  CL (2)  (x)  x  dx 

Zn 


rr  n=Q        n       r    dr   •-'    r         J 

It  can  be  shown  that 

(2 
2n 


^  2 

x    -     r 


/7T7ci2»(X)  = 


1 


-  (2n  +  1)U 


2n  +  2 


4 
(x)] 


=.[(2n+3)U_      (x) 

J^1  2n 


(ID 


(12) 


and  therefore 


N 


g(r)=+£      a     [(2n+3)I      (r)-(2n+l)I      ,    .   (r)]  , 

n=0     n  n  n  +  1 


(13) 


where 


,       .    U.     (x)xdx 
T     i    i-"1      1 d    pi       2n 

nKT}      2n     r     dr    J  r     f~     ~T   I      2        T 

yi   -  x      7x       -  r 


and  thus  I     (r)  =  0  and  I.  (r)  =   -  1. 
o  I 


(14) 
(15) 


Siibstituting  the  recurrence  relation 

U  (x)=  2(2x2  -  1)  IT        ■        (x)  -  IT        (x) 

2n  +  4  2n  +  2  2n 


=  +2U  (x)-U         (x)  +  4(x     -1)U  (x) 

2n  +  2  2n  2n  +  2 


(16) 


into  Eq.(14),we  obtain 


2      1       d     ,.1 


1  2 

yi-x  u 


n 


+  _=2I  -I    +^-L-iLJ 

+  2  n+1       n       tt      r   dr    ,J 


2n  +  2 


(x)  x  dx 


n — 


(17) 


The  integral  of  the  last  term  has  a  closed  form  solution  (see  Table  1) 
so  that  Eq.    (17)  becomes 


J„+  2(r)=   2  I„+l<r»  "  yr»  -  <2n+  3)Pn+  l(2r2  -  ,  ,,  (18) 

2  5 

where    P  ,(2r     -1)  is  also  generated  by  a  recurrence  relation 

n  +  1 

given  by 


^.l^-1'^^-1)^2-1' 


P  ,(2r2  -   1).  (19) 


n  +  1       n  -   1 

Using  Eqs.    (15),    (18),    and  (19),    the  desired  function,    g(r),    can 
readily  be  obtained  from  Eq.    (13). 

This  approach  was  readily  programmed  for  a  computer  solution. 

The  integral  of  Eq.    (10)  was  evaluated  using  the    trapezoidal  rule  and  a 

40  point  Gaussian  quadrature.      The  change  in  the  values  of  the  coefficients 

using  the  trapezoidal  rule  was  completely  negligible.      A  listing  of  the 

Fortran  computer  program  using  the  trapezoidal  rule  is  given  in  Appendix  II. 

The  series  of  coefficients  so  calculated  can  be  terminated  by  applying 

the    F-distribution  significance  test     to  each  coefficient  in  turn.     We 

start  by  assuming   that  a     is  always  significant.     We  then  calculate 

the  mean  square  deviation,    6,    between  the  input  curve    f(x.)  and  the 

l 

curve  defined  by  Eq.    (2)  using  the  coefficients  known  to  be  significant. 
This  calculation,    when  repeated  using  one  additional  coefficient  to  be 
tested,    yields     6  '.      The  quantity  F*  =  (6   -  6  ')  x  (number  of  points  in  the 
input  curve )/5  '     is  then  compared  with  the  value  of    F     chosen  from 


10 


Table  V  of  Ref,    6.     If  F*  >  F    then  the  additional  coefficient  is  considered 
significant  and  the  procedure  is  repeated  for  the  following  coefficient.    When 
a  non- significant  coefficient  is  found,    all  subsequent  coefficients  are 
then  set  equal  to  zero.      This  is  done  on  the  physical  basis  that  fine 
structure  is  not  expected. 

In  treating  the  input  curves  we  allowed  for  asymmetries  by 
following  the  approach  of  Freeman  and  Katz,      and  thus  did  not  force 
the  transformed  curves  to  be  symmetric.      This  was  done,    in  spite  of 
the  fact  that  no  large  asymmetries  were  expected,    in   order  to  obtain 
a  check  on  the  symmetry  of  the  discharge.      In  this  approach,    the 
symmetric  part  of  f(x)  is  contained  in  the  even  function 


f  (x)=  M+*)  +  f(-x)  (20) 

g  z 

and  the  asymmetric  part  of    f(x)  is  described  by  the  even  function 

f  (+  x)  -  f  (-x)  (21) 

u  2x 

After  the  Abel  inversion  of  these  functions  to  obtain  g   (r)  and  g    (r), 

g 

the  radial  distribution  is  constructed  by  using 


g  (+  r)=  g   (r)  +  rg   (r)  (22) 

8 


and 


g(-r)  =  g   (r)-rg   (r).  (23) 

6 


11 


CONCLUSION 
The  use  and  advantages  of  orthogonal  polynomial  expansions  in 
solving  the  Abel  integral  transform  have  been  briefly  discussed  and  a 

new  solution  in  terms  of  Gegenbauer  (a  =  2)  polynomials,    which  utilizes 

2 

a  weighting  function  (1-x    )  in  the  calculation  of  the  expansion  coefficients, 

has  been  presented.     In  the  presence  of  significant  distortions  near  the 

boundary  of  the  luminous  region,    a  condition  which  is  particularly  true 

2 

for  projected  profiles  approximating  f(x)  =       1-x   ,    or  g(r)  approximately 

constant,    the  use  of  the  Gegenbauer  (a=2)  polynomial  expansion  will 
reduce  the  contribution  of  such  distortions  in  the  calculation  of  the 
expansion  coefficients,    and  in  principle  allow  a  more  accurate  solution 
of  the  Abel  inversion  integral. 
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APPENDIX  I 
We  present  here  a  proof  of  the  statement  that  the  "impact 
parameter"  of  a  light  ray  is  unchanged  by  passage  thru  the  wall 
of  a  uniform  glass  tube.     Referring  to  figure  1,    we  write  the  law  of 
sines  for  triangle  ABC  as  sin  0     =   sin  ( 1 80-0    )/( 1  +  t  ),    or 


(1  +  t)  sin  92  =   sin  0    .  (A-l) 


The  initial  and  final  impact  parameters  are  given  by  the  relations 


p.=  (l+f)sin9.  ,  and  p    =   s  in  9  „ .  (A-Z) 

*i  1     .         rf  4 


The  angles  involved  are  related  by  Snell' s  law  as  follows: 

sin  9    /  sin  9     =  sin  8    /  sin  9     =  r|.  (A-3) 

Here    r|     is  the  index  of  refraction  of  the  glass  tube  relative  to  that  of 
air  and  the  indices  of  refraction  both  inside  and  outside  of  the  tube  are 
assumed  identical.      Combining  these  relations,    we  obtain 
p    =   sin  9 -.  -   sin  0      (sin  9    /sin  9    )  =   sin  9      (1  +  t)  =  p..      Thus  the  impact 
parameter  inside  the  tube  is  identical  with  that  of  the  incident  ray.      This 
result  is  simply  another  expression  of  the  conservation  of  angular 
momentum.      There  is  consequently  no  distortion  due  to    refraction  for 
the  case  of  a  cylindrically  symmetric  light  source,    coaxial  with  a 
cylindrical  glass  tube  of  uniform  thickness. 
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INCIDENT 
LIGHT   RAY 


Figure  1.      Displacement  of  an  incident  beam  of  light  in  passing 
through  an  optically  perfect  glass  tube.      The  "impact 
parameter"   is  not  changed. 
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IB 


APPENDIX  II 


PROGRAM  ABELTAPE 

DIMENSION  HFD(  5 ) ♦ X ( 500 ) ♦ Y( 500 )  ,RP ( 100 ) ,RM(  100 ) »XP ( 250 ) »YP ( 250  )  , 
1YM(?50) »XF(501) »FX(501) •URK501 ) »UR2(501 ) *YKN(25l ) »AN(?0) »BN(20) » 
2ANN(20) »BNN(20) ,RA ( 20 ) »RB ( 20 ) »RX ( 5  01 ) » CS ( 501 )  »BS ( 501 ) »GP ( 100) » 
3GMQ00)  » KIND (3,2  )>GR(201)  »XJP(251)  »XJM(251)  »XX(251  )  ,HEE(10) 
TYPE  INTEGER  DATE 
COMMON  NPOINT 
DATA(  IRUNO  =  0)  »( ISETO=0) 

DATA (ESP  I  =  5. 092958)* ( FSP= 1. 273  2  395 )  »( (KIND(ltl) »I=1»2)=8H    TCHEB, 

15HICHEF) ,  (  (KIND(2»I ) ,I=1»2)=8H     LEGE.4HNDRE)  »((KIND(3»I)»I=1»2)= 

?8H     GEGE.5HNBAUR) 

112  FORMAT(*0SET=*,I2»5X,*RUN=*»I2»*   NOT  FOUND  ON  INPUT  TAPE*) 

READ(60,505)HEE 
505  FORMATM0A8) 
C 

C  GENERATE  RHO  ARRAY 

C 

RZ  =  0. 
RP<1)=.01 
RM(1)=-.01 
DO  106  I=2»100 
RP(  I  )=RP(  I-D  +  .01 
106  RM( I )=-RP(  I  ) 

DATE=IDATE(XDUMMY) 
C 

C  GENERATE  GR  ARRAY 

C 

DO  500  i=i,ioo 

IM=101-I 
IP=101+I 
GR( I )=RM( IM) 

500  GR( IP)=RP(  I  ) 
GR( 101 )=0. 

C 

C  READ  INPUT  PARAME TERS,DATA »  AND  HEADING. 

C 

99  READ (60,1 11 ) I  SET , I  RUN ,NN ♦ NX »HED  ♦  I NFLAG 
111  F0RMATUI5  ,5A8,I1) 
IF(EOF, 60)98, 501 

501  IF( INFLAG.EQ.O)   GO  TO  503 
READ(60,153) ( X ( I ) ♦ Y( I ) , 1=1 ,NX ) 

153  FORMAT(8(F5.0»F5.1  )  ) 

GO  TO  9999 
5  03  READ (60,151) ( X ( I ) » 1  =  1 »NX ) * ( Y(I ) ♦ 1= 1  »NX ) 

151  FORMAT(8E10.3) 
999  9  READ(60,152)NPOINT,NP,FVALUE,THICK»IPPP 

152  FORMAT(2I5»2F10.5»  II  ) 
C 

C  GENERATE  ARRAYS  FOR  YM  AND  YP ♦  THE  CORRESPONDING  Y 

C  VALUES  FOR  -X  AND  +X  RESPECTIVELY  AND  ASSUME  THE 

C  VALUES  AT  THE  END  POINTS  ARE  AVERAGFS  OF  THE  MEASURED  Y. 

C  NH=NX/2  FOR  EVEN  NX. 

C  NH=(NX-l)/2  FOR  ODD  NX. 

C  NT  IS  THE  TOTAL  NO.  OF  POINTS. 

16 


12 


14 


13 


15 


16 


FIRST  FIND  THE  AVERAGE  DELTA  X# 


6  DXA=2./FL0ATF(NX-1 ) 


FIND  N=  NX   MODULO   2. 

N=XMODF(NX*2) 
IF(N.EQ.O) 12»13 

N=0  IMPLIES  NX  IS  EVFN. 

NH=NX/2 

DXH=DXA/2. 

NHPO=NH+l 

DO  14  I=1»NH 

NHP=NH+I 

YP( I )=Y(NHP) 

XP(  I  )=DXA*FLOATF( 1-1 )+DXH 

NHM=NHPO-I 

YM( I )=Y(NHM) 

YZ=.5*( Y(NH)+Y(NHPO) ) 

GO  TO  16 

N  NOT  0  IMPLIES  NX  IS  ODD, 

NH=(NX-1 )/2 

NHPO=NH+l 

DO  15  1=1 »NH 

NHP=NHPO+I 

YP (  I) =Y(NHP) 

NHM=NHPO-I 

XP(  I  )=DXA*FLOATF(  I  ) 

YM( I )=Y(NHM) 

YZ=Y(NHPO) 

DXH=DXA 

XP(NH)=1. 

YA=.5*( YP( NH)+YM(NH)  ) 

YP(NH)=YA 

YM(NH)=YA 


APPLY  BACKGROUND  NOISE  CORRECTION  TO  THE  DATA, 


C  =  0. 

C=YA/SQRT(TPS-1.) 


17 


TP=THICK+1. 

TPS=TP*TP 

IF (THICK. EQ.ru  ) 

IF(THICK.NE.O. ) 

YZ=YZ-C*THICK 

DO  17  I=1,MH 

XSQ=XP(  I  )*XP(  I ) 

DXP=SQRTF( TPS-XSQ)-SQRTF( l.-XSQ) 

CD=C*DXP 

YP(  I  ) =YP(  I  )-CD 

YM(  I  )=YM(  I  )-CD 


17 


19 
21 
22 
18 


20 


23 


TRANSFORM  THE  FUNCTION  SO  THAT  IT  IS  ZERO  AT  THE 
BOUNDARY  AND  IT  IS  NORMALIZED. 

YZ=YZ-YA 

YMAXIN=YZ 

DO  18  1  =  1  ,NH 

YP(  I  )=YP(  I  )-YA 

YVK  I  )=YM( I )-YA 

IF(YP(I).GT,YMAXIN)19,21 

YMAXIN=YP( I ) 

IF(YM( I ) .GT. YMAXIN )22»18 

YMAXIN=YM( I ) 

CONTINUE 

YZ=YZ/YMAXIN 

DO  20  I=1»NH 

YP(  I  )=YP(  I )  /YMAXIN 

YM(  I  ) =YM(  I ) /YMAXIN 


24 


SET  UP  XF  AND  FX  ARRAYS  USED  RY 
SIGNIFICANT  COFFFICIENT  TEST. 


HE 


NHPO=NH+l 

NT=NH+NHPO 

DO  23  I=1»NH 

NHM=NHPO-I 

FX( I )=YM(NHM) 

XF(  I  )=-XP(NHM) 

NHP=NHPO+I 

FX (NHP)=YP  (  I  ) 

XE(NHP)=XP(  I  ) 

FX(NHPO) =YZ 

XF(NHPO)=0. 

CALCULATE  SYMMETRIC  AND  ASYMMETRIC  PORTIONS  OF  FUNCTION. 

NHPO=NH+l 

XJP( 1)=YZ 

XJM( 1)=0. 

DO  24  I=2»NHPO 

XJP( I )=.5*(YP( I~1)+YM( 1-1 ) ) 

XJM(  I  )  =  (YP(  I.-1)-YM'(I-1  )  >/(2.*XP(  I'-l)  ) 

WRITE  OUT  FIRST  PAGE  OF  OUTPUT. 


WRITE(61»121)HEE»HED»DATE»ISET» I  RUN ,NN ,NX * F VALUE , THI CK » NP , K IND ( NPC 
1INT»1) •KIND(NPOINT»2) 

121  FORMAT(  1H1  ♦ 1 0A8 » 5A 8 » 3X ♦ A8  //6H  SET  =  I  3  ,4X  *  4HRUN=  I  3  »V+) 
l»*TOTAL  NO  POINTS  =  *» I4»4X ,*NO  POINTS  USED=* *  I  4 ♦ 4X • 8HF  VALUE  =  F6.2  ♦'; 
2X,6HTHICK=F8.3,4X,3HFORI3»lX,2A8,6H  TERMS//) 

WRITE (61  ,122  ) 

122  EORMAT(*  THE  LISTING  OF  THF  IN^UT  DATA  FOLLOWS*//) 
WRITF  (61  ,123) 

123  FORMAT( 5 ( 1  OX ♦ 1 HX »8X , 1HY ♦ 2 X  )  ) 
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WRITE(61 tl 

124  F0RMAT(5(F 
WRITE(61,1 

1INT,1 ) »KIN 
WRITE (61 »1 

125  FORMAT(»  T 
WRITE  (61  ».l 
WRITE (61 , 1 
GENERATION 
DXH=DXH/2. 
XX(1 )=n. 
DO  91  I=2» 

91  XX(  I )=XP(  I 
GO  TO  (100 

100  DO  101  1=1 
UR2( I )=1. 
URK  I  )=XX( 
IF( I.FQ.ll 

19  8  AN(1)=XJP( 
BN( 1 ) =  XJM( 
GO  TO  101 

199  AN( 1 )=AN(1 
BN(1 )=BN  (1 

101  CONTINUE 
AN( 1 ) =AN  (1 
BN( 1 ) =BN  (  1 
DO  102  J=2 
DO  103  1=1 
TX=XX(  I  >  +  X 
UR2( I )=TX* 
UR1 ( I )=TX* 
IF(  I.EQ.l) 

104  URR=UR2( I ) 
AN( J)=XJP( 
BN( J)=XJM( 
GO  TO  103 

105  URR=UR2( I ) 
AN( J)=AN( J 
BN( J)=BN( J 

103  CONTINUE 

AN( J)=AN (J 
BN ( J) =BN (J 

102  CONTINUE 
GO  TO  n 

200  DO  201  1=1 
UR2( I  )  =  1. 
YKN( I)=2.* 
UR1 ( I )=YKN 
IF( I.EQ.l) 

298  URR=DXH*XX 
AN( 1 )=XJP( 
BN( 1 )=XJM( 


24)  (X(  I )  »Y(  I  )  ,1=1, NN) 
13.4,F9.4)  ) 

21)HEE,HFD,DATE,ISFT, I  RUN ,NN ,NX , FVALUE , THI CK ,NP , K IND ( NPO 
D(NPOINT,2) 
25) 

HE  LISTING  OF  THF  NORMALIZED  INPUT  DATA  FOLLOWS*//) 
23) 

24)  (XF( I  )  ,FX(  I )  ♦  1  =  1, NT) 
OF  THE  COEFFICIENTS  FOR  THE  DESIRED  POLYNOMIAL 


NHPO 

-1) 

,200,300) ,NPOINT 

NPOINT=l, 
,NHPO 


TCHEBICHEF  CASE 


I  )+XX(  I ) 
198,199 
1  )*DXH 
1 )*DXH 

)+XJP( I )*DXA 
)+XJM( I )*DXA 

)*FSP 

)*FSP 

,NP 

,NHPO 

X(  I  ) 

UR1  (  I  )-UR2( I ) 

UR2(  I  )-URl ( I ) 

1 04,105 

*DXH 

I  )*URR 

I  )*URR 

*DXA 

)+XJP( I )*URR 

)+XJM(  I  )*URR 

)*FSP 
)*FSP 


,NHPO 

XX(  I  )*XX(  I  )-l. 

(  I  ) 

298,299 

(  1) 

1 )*URR 

1)*URR 


NPOINT=2,  LEGFMDRF  CASE 


19 


299 

201 

234 

235 
231 

244 

245 


249 

204 


205 


203 

255 
254 


20? 


30n 


398 


399 


301 


GO  TO 
URR  =  D 
AN(1  ) 
BN(  1) 
CONTI 
AN(  1  ) 
BN(  1  ) 
DO  20 
IF(J. 
XN  =  J- 
GO  TO 
XN  =  J- 
XN1  =  J 
XN2  =  X 
DO  20 
IF(  J. 
UR  =  UR 
GO  TO 
UR=(X 
UR2(  I 
URK  I 
IF(  I. 
URR=U 
AN  (J) 
BN(  J) 
GO  TO 
URR=U 
AN(J) 
BN(  J) 
CONTI 
IF(  J. 
XN2  =  X 
TX  =  XN 
AN(  J) 
RN(  J) 
CONTI 
GO  TO 

DO  30 
YKN(  I 
UR2  (  I 
UR1  (  I 
IF(  I. 
URR  =  Y 
AN(  1  ) 
BN(1) 
GO  TO 
URR  =  Y 
AN(1) 
BN(  1  ) 
CONTI 
TX  =  ES 
AN  (  1  ) 
BN(  1  ) 


201 
XA*XX 
=  AN  (1 
=  BN  (1 
NUE 
=  2.*A 
=  2.*B 

2  J  =  2 
EQ.2) 
1 

231 
2 

-1 
N  +  XN  + 

3  1  =  1 
EQ.2) 
KI) 

249 
N2*YK 
)=UR1 
)=UR 
FQ.]  ) 
R*DXH 
=  XJP( 
=  XJM( 

203 
R*DXA 
=  AN  (J 
=  BN  (J 
NUE 
E0.2) 
N2  +  2. 
2*2. 
=  AN(  J 
=  BN  (J 
NUE 

n 

i  i  =  i 
)  =  i.- 
)  =  i. 

)=4.* 

EQ.l  ) 
KN(  1) 
=  XJP( 
=  XJM( 

301 
KN(  I  ) 
=  AN(] 
=  BN  (1 
NUE 
PI/3. 
=  AN  (1 
=  BN(  1 


(  I  ) 

)+XJP{ I )*URR 

)+XJM( I )*URR 

N(  1  ) 
N(  1) 
,NP 
234,235 


1. 

♦  NHPO 

244*245 


N(  I )*UR1  (  I  ) 
(I  ) 

204,205 
*XX(  I  ) 
I )*URR 
I  )*URR 

*XX(  I  ) 

)+XJP( I ) *URR 
)+XJM( I )*URR 

254,255 


)*TX 
)*TX 


,NHPO 

XX(  I  )*XX(  I  ) 

XX(  I  ) 
398,39^ 
*DXH 
1 )*URR 
1  )*URR 

*DXA 

)+XJP( I ) *URR 

)+X JM( I )*URR 


•XN*UR2 ( I ) )/XNl 


NPOINT=3,  GEGENBAUR  CASE 


)*TX 
)*TX 


20 


304 


305 


303 


302 


11 


31 


33 
34 

32 
36 
35 

38 


41 

39 
4? 

51 


52 


XN  =  0. 
DO  30 
XN  =  XN 
XN1=X 
XN2  =  X 
XN3  =  X 
XN4  =  X 
TXRAT 
DO  30 
UR2(  I 
UR1  (  I 
UR=UR 
IF(  I. 
URR=U 
AN(  J) 
BN(  J) 
GO  TO 
URR=U 
AN(  J) 
BN(J) 
CONTI 
TX  =  ES 
AN(  J) 
BN(  J) 
CONTI 


2  J 
+  1. 
N  +  l 
N1  + 
N2  + 
N3  + 
=  XN 

3  I 
)  =  ( 
)  =  ( 
2(  I 
EQ. 
R*D 
=  XJ 
=  XJ 

30 
R*D 
=  AN 
=  BN 
NUE 
PI/ 
=  AN 
=  BN 
NUE 


■•?»NP 


1. 

1. 

1  . 

2* 

=  1 

2. 

2. 

)* 

1) 

XH 

P( 

M( 

3 

XA 

(J 

(J 


XN4 

♦  NHPO 

*XN2*  XX ( I )*UR1 ( I )-XN3*UR2 ( I ) ) /XN1 

*XN3*  XX ( I )*UR2( I )-XN4*URK I ) )/XN2 

YKN(  I  ) 
304,305 

I )*URR 
I )*URR 


)+XJP( I )*URR 
)+XJM( I ) *URR 


TXRAT 
( J)*TX 
( J)*TX 


TEST  FOR  SIGNIFICANT  COEFFICIENTS. 


DO  31  I=1»NP 

ANN( I )=AN( I ) 

BNN( I )=BN( I ) 

RA(  I  )=0. 

RB( I )=0. 

RA(1)=FVALUE+1. 

GO  TO  (32»33»32) *NPOINT 

DO  34  K=l »NT 

RX(K)=1. 

GO  TO  35 

DO  36  K=1>NT 

RX(K)=SORTF( l.-XF(K )*XF(K ) ) 

DO  37  1=1, NP 

IF(  I.EQ.l) 38,39 

DO  41  J=l ,NT 

UR1 ( J)=l. 

UR2 ( J)=l. 

BS( J)=AN (1 ) 

GO  TO  26 

IF( I.EQ.2) 43,44 

DO  45  J=1 ,NT 

GO  TO  ( 51  ,52,5?)  ,NPOINT 

UR2(J)=4.*XF(J)*XF{J)-1. 

UR1 ( J)=XF( J)+XF( J) 

GO  TO  54 

Z=2.*XF( J)*XF( J)-l. 
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53 

54 
45 

44 

61 


62 


63 


64 
47 
46 


48 


26 


68 


37 


65 


66 


UR2( J)=Z 

GO  TO  54 
UR2( J)=12* 
UR1 ( J)=4.* 
CS(  J)=8M  J 
BS( J)=BS( J 
GO  TO  46 
DO  47  J=l, 
GO  TO  (61» 
TXX=XF( J)+ 
URK  J)=TXX 
UR2(J)=TXX 
GO  TO  64 
XM2  =  I 
XMN=I-1 
XM1=XM2+XM 
Z=2.*XF( J) 
TXX=(Z*UR2 
URK  J)=UR2 
UR2( J)=TXX 
GO  TO  64 
XM2=I+I-3 
XM3=XM2+1. 
XM4=XM3+1 . 
XM5=XM4+1 • 
URK  J)  =  (2. 
UR2( J)=(2. 
CS( J)=BS( J 
BS( J)=BS( J 
AS  =  0. 
DS  =  0. 

DO  48  K=l» 
ASR=FX(K)- 
DSR=FX(K)- 
AS=AS+ASR* 
DS=DS+DSR* 
RA(  I  )=ABSF 
AS  =  0. 

Ds  =  n. 

DO  68  K=l ♦ 

CS(K)=BS(K 

BS(K)=BS(K 

ASR=FX(K)- 

DSR=FX(K)- 

AS=AS+ASR* 

DS=DS+DSR* 

RB( I )=ARSF 

CONTINUE 

IPP  =  1 

DO  42  1=1 » 

IF(RA(  I  )  .L 

DO  66  J=I, 

ANN( J)=n. 

BNN( J)=n. 


XF( J)*XF( J)-2. 

XF(  J) 

) 

)+AN( I )*UR2( J)*RX (J) 

NT 

62»63) »N<=OINT 
XF(  J) 

*UR2(J)-UR1( J) 
*UR1 ( J)-UR2( J) 


N 

*XF( J)-l. 

( J)*XM1-XMN*UR1 ( J) ) /XM2 

(J) 


*XM3*XF( J)*UR2( J)-XM4*UR1 (J) )/XM2 

*XM4*X^( J)*UR1 ( J)-XM5*UR2( J) ) /XM3 

) 

)+AN( I )*UR2( J)*RX ( J) 


NT 
CS(K) 

BS(K) 

ASR 

DSR 

( 100.*(AS-DS)/DS) 


NT 

) 

)+BN( I )*UR2(K)*XF(K)*RX( K) 

CS(K) 

BS(K) 

ASR 

D.SR 

(10n.*(AS-DS) /DS) 


NP 

E.FVALUE )65*69 

NP 
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GO  TO  85 
69  GO  TO  (84,42)  » I PP 

84  IF(RB(I ).LE.FVALUE)71,42 

71  DO  72  J=I,NP 

72  BNN(J)=0. 
IPP  =  2 

42  CONTINUF 

85  IPP=1 
C 

C  GENERATE  F(RHO)  VERSUS  RHO. 

C 

67  CALL  FGRHO(RZ,ANN,l,GZ,NP  ) 

DO  73  1=1,100 

CALL  FGRHO(RP( I )  » ANN  *  1 ,GP ( I )  »NP ) 

CALL  FGRHO(RM( I ) » PNN , 2 »GM ( I ) » NP ) 

73  CONTINUF 
C 

C  WRITE  OUT  SFCOND  PAGF  OF  OUTPUT. 

C 

WRITE(61»121)HEE»HFD,DATE,ISFT, I RUN »NN »NX ,F VALUE, THICK »NP ,K IND ( NPO 
1INT,1 ) »KIND(NP0INT,2 ) 

WRITE(61,126) 

126  F0RMAT(8X»*SIG  COF FF*  ,21 X  ,*ALL  COEFF* , 24X ,*RAT I0S*/4X ,*A ( I ) * ,9X ,*B 
1(  I  )*,13X,*A(  I  )*»9X»*B(  I)*»13X»*RA(  I  )*»10X,*RB(  I  )*/'/) 

WRITE (61, 127)  (ANN(  I  )  ,RNN(  I )  , AN (  I )  , BN ( I  )  , RA (  I  )  , RR (  I  ) ♦ I  ,  I  =  1  ,NP ) 

127  F0RMAT(F10.6,F13.6,F17.6,F13.6,E19.4 ,E15 »4, 6X » I  2 ) 
C 

C  NORMALIZE  F(RHO)  AND  SFT  TO  ZERO  IF  IT  IS  NEGATIVE. 

C 

81  DO  74  1=1,100 

IM=101-I 

IP=101+I 

FX(  I)=GP( IM)-GM( IM) 

74  FX( IP)=GP( I )+GM( I ) 
FXdOl  )=GZ 
FMXOUT=FX( 1 ) 

DO  75  1=2,201 

IF(FX( I ) .GT.FMXOUT) 76,75 

76  FMXOUT=FX( I ) 

75  CONTINUE 

83  DO  77  1=1 ,20] 

FX( I )=FX( I ) /FMXOUT 

IF(FX( I ) .LT.0.)78,77 
78  FX( I)=0. 

77  CONTINUE 
C 

C  WRITE  OUT  THIRD  PAGE  OF  OUTPUT. 

C 

WRITE(61,121 )HEE,HED,DATE,lSET,IRUN,NN,NX,FVALUE,THICK,NP,KIND(NPO 

1INT,1 ) ,KIND(NPOINT,2  ) 
WRITE(61,128) 

128  FORMAT(5(9X,3HRHO,4X,6HF(RHO) ) //) 
WRITE (61, 124) (GR(I ) ,FX( I) ,1=1,201) 

WRITE (61 , 1 29 )YMAX IN, FMXOUT 
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129  FORMAT(*OTHE  NORMALIZATION  FACTORS  ARE*//*  INPUT/* ,E 1 1 .4 » 10X »*OUTP 
1UT/*»E11.4) 

IF( IPPP.EQ.l ) 9999,99 
98  WRITE(61»131  ) 
131  FORMAT(*  END  OF  THIS  RUN*) 

WRITE(61 ,506)HEE,DATE 
506  FORMATdHl  ♦ 1 0A8 ♦ 20X , A8 ) 

CALL  EXIT 

END 

SUBROUTINE  FGRHO ( X ♦ AN ♦ J »FOX »NP ) 

DIMENSION  AN(20) 

COMMON  NPOINT 

DATA(TP=6. 2831853) 

Y=2.*X*X-1. 

XX  =  3. 

GO  TO  (  11»12*13)  »NPOINT 

11  P2=l. 
P1=Y 

S=AN(   1)+3.*AN(   2)*Y 

XN  =  0. 

DO  1  I=3»NP 

XX=XX+2. 

XN=XN+1. 

PN=( (2.*XN+1. )*Y*P1-XN*P2 )/ (XN+1. ) 

S=S+XX*AN(   I)*PN 

P2  =  P1 

P1  =  PN 

1  CONTINUE 

GO  TO  (2*3  )  »J 

2  FOX=.5*S 
RFTURN 

3  FOX=-.^*S*X 
RETURN 

12  XSQ=X*X 
ZSQ=1.-XSQ 
TZSQ=ZSQ+ZSO 
Z=SQRTF(ZSQ) 
P2  =  l. 
P1=3.-4.*XSQ 

S=AN(   1)-AN(   2)*P1 

S  I  =  1  . 

DO  4  I=3»NP 

PN=TZ9Q*P1-P2 

P2=P1 

P1  =  PN 

S=S+SI*AN(    I)*PN 

4  SI=-SI 
DEN=1,5707963*Z 
GO  TO  (5*6)  »J 

5  FOX=S/DEN 
RETURN 

6  FOX=-S*X/DEN 
RETURN 

13  P?=1. 


24 


P1  =  Y 
T2=-l. 

T1=1.-6,*X*X 

S  =  AN(       1)+AN(       2)*(5.*T2-3.*T1  ) 

DO    8    I=3»NP 

XN=I-1 

XNN=XN-1. 

XB=XN+XN+1. 

XA=XB+2. 

PN=( (XN+XNN)*Y*P1-XNN*P2) /XN 

TN=T1+T1-XB*PN-T2 

I )*(XA*T1-XB*TN) 


S=S+AN( 

I  ) 

*( 

T2  =  T1 

T1  =  TN 

P2  =  P1 

P1=PN 

GO  TO  (7 

.9) 

♦  J 

FOX=    S" 

*X 

RETURN 

FOX  =  S 

RETURN 

END 
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